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Abstract 

Using the tools of semiconvex duality and max-plus algebra, this work 
derives a new fundamental solution for the matrix differential Riccati 
equation (DRE) with time-varying coefficients. Such a fundamental solu- 
tion, is the counterpart of the state transition matrix in linear time- varying 
differential equations, and can solve the DRE analytically starting with 
any initial condition. By parametrizing the exit cost of the underlying 
optimal control problem using an additional variable, a bivariate DRE is 
derived. Any particular solution of such bivariate time-varying DRE, can 
generate the fundamental solution, and hence the general solution, ana- 
lytically. The fundamental solution is equivalently represented by three 
matrices, and the solution for any initial condition is obtained by a few 
matrix operations on the initial condition. It covers the special case of 
time invariant DRE, and derives the kernel matching conditions for trans- 
forming the DRE into the semiconvex dual DRE. As a special case, this 
dual DRE can be made linear, and is thus solvable analytically. Using 
this, the paper rederives the analytical solutions previously obtained by 
Leipnik [4] and Rusnak [6]. It also suggests a modification to the expo- 
nentially fast doubling algorithm described in [1], used to solve the time 
invariant DRE , and makes it more stable and accurate numerically for the 
propagation at small time step. This work is inspired from the previous 
work by McEneaney and Fleming [1],[2]. 
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1 Introduction 



In this paper, we consider the matrix differential Riccati equations (DRE's) of 
the form 

= A{typ{t) + p{t)A{t) + c{t) + p{tymp{t) 

p{T) ^pt and t<T ^ ^ 

where t € R and A{t) is square and p(t), C(t), E(i) are square and symmetric 
nxn matrices. Note that (f ) can easily by converted to an initial value problem 
with forward time propagation, but we use this approach because it simplifies 
the notation when framed as an optimal control problem. 

DREs are widely used in system and control theory, specifically in optimal 
control, filtering and estimation. Many numerical algorithms have been pro- 
posed in past for solving time- varying DREs. These include carefully redesigned 
conventional Runge-Kutta methods and other explicit linear Multi-step meth- 
ods for ODEs, as well as nonlinear implicit by Choi and Laub [11], Dieci [12] and 
many others. Although these methods benefit greatly from past development in 
general purpose computer programs for solving ODEs, they can become rather 
complex in code structure and interface. Implicit methods, which are more 
preferred to the explicit ones for solving stiff problems, also suffer from imple- 
mentation and computational complexity. Also these methods have to be rerun 
to solve for each initial condition, making sensitivity analysis rather difficult. 

There have also been unconventional methods for the DREs arising in opti- 
mal control, e.g. [5], [6], [9], [13], [14] and [15]. These cover various analytical 
solutions, and doubling algorithms for time-invariant problem. But it is known 
that these are not suited for time-varying DREs. One method which can be 
used is the analytical solution by Davison and Maki [17]. It solves the following 
system. 

Pt,=V{h)U-\h) (2) 

with the solution obtained as Pt^ = V{t2)U^^{t2). Thus the method does work 
for the time- varying systems. But as ^2 — grows, columns of U{t2) become 
more and more linearly dependent, making the problem ill-conditioned, hence 
it can be used only for small time propagation. Thus until now, there has been 
no fundamental solution available for the time-varying DREs, which is useful 
for the long time horizon propagation, and for the infinite horizon as a special 
case. This work attempts to fill this gap. 

Similar to the method by Sorine and Winternitz [9], this paper provides a 
way to construct a general solution from the particular solutions. But instead of 
using five particular solutions starting at rather special initial conditions as in 
[9], this paper uses just one particular solution of the bivariate DRE, starting at 
any initial condition to construct the general solution of the time varying DRE. 

Now we shall discuss an overview of the forthcoming development. Re- 
cently McEneaney [1], proposed a new fundamental solution for solving the 
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time-invariant DRE, and a doubling algorithm. They are based on the tools 
of max-plus algebra and semiconvex duality. Similarly [2] and [8] introduced 
the concept of the max-plus fundamental solution, for the time invariant sys- 
tems with nonlinear dynamics. This paper extends both of these ideas to apply 
to a time-varying system with linear dynamics and quadratic payoff. It finds 
surprisingly simple formulae for the fundamental solution, which turns out to 
be bivariate quadratic and can be represented by three n dimensional square 
matrices. It can also be viewed as the max-plus kernel, which can operate on 
any initial condition po to obtain pt analytically. 

The fundamental solution is obtained from the time evolution of a bi-variate 
quadratic terminal cost function under optimal control, where the second vari- 
able is used to parametrize the terminal cost, the other variable being the state. 
Note that this requires us to solve a bi-variate DRE (11), hence evolve three 
parameters instead of just one in case of the DRE (1). The fundamental so- 
lution itself is invariant with respect to the terminal bi-variate cost function, 
but depends only on evolution time interval. Thus any particular solution to a 
bivariate DRE can be converted into the fundamental solution, and thus into a 
general solution, analytically. This makes sensitivity/ perturbation analysis for 
such initial value problem much easier. 

As a special case, the fundamental solution to the time-invariant problem 
and a new doubling algorithm is derived. Though inspired from the doubling 
algorithm in [1], the new algorithm is more direct and simpler, since it does away 
with the kernel propagation in semiconvex dual space. Instead, it does so in the 
primal space of quadratics. One numerical issue with such fundamental solution 
as well as the one described in [1], is that the kernel parameters blow up as the 
time step gets smaller. A modification to the formula to avoid such blowup is 
suggested, which maintains better solution accuracy for the propagation at a 
small time step. 

Building further on the idea of semiconvex duality, it shall be proved that 
the semiconvex dual of the solution of DRE (1) satisfies another DRE whose 
coefficients can be found analytically. One can also choose these coefficients, 
and find an appropriate duality kernel for transformation, as long as certain 
compatibility/matching conditions are satisfied. Such conditions take the form 
of coupled Riccati equations or similarity transformation on the primal and dual 
Hamiltonian matrices. As a special case, one can make the dual DRE linear, by 
choosing the quadratic term coefficient zero, and solve it analytically. Using this 
method, the analytical solutions obtained earlier by Leipnik [4] and Rusnak [6] 
can be easily derived, demonstrating the versatility and power of this approach. 

We shall also use the fundamental solution developed here to solve stiff DREs 
with known analytical solutions, and benchmark them for accuracy, numerical 
stability and speed. These algorithms, being analytical and stable, are very 
useful for solving the stiff DRE for long time horizon propagation, unlike the 
Davison-Maki method. 

Finally even though the optimal control problem considered here does make 
a number of assumption to ensure existence of the value function, and avoid 
singularity, resulting fundamental solution is valid for a much wider class of 



problems, since it is purely algebraic. As a special case, it is observed that even 
in case of unstable DREs, which exhibit finite time blowup, above fundamental 
solution can propagate beyond the singularity, e.g. when applied to special 
case p = 1+ p^, whose solution is tan(p), it correctly propagates beyond the 
singularities at p = (2n + 1) f , n G W, where W is the set of whole numbers. 
But this aspect shall not be covered here in order to contain the scope of the 
current paper. 



2 Optimal control problem 

We shall obtain the fundamental solution for DRE (1) through an associated op- 
timal control problem. To ensure existence and regularity of the value function, 
we make following assumptions throughout this section. 

Since DREs exhibit finite time blowup, we assume that there exists 

a solution of DRE (1) ,with terminal condition Pt for t € (f, T] 

with t<T. We may have T = -oo. We assume that S(t) ^ S 

{f,T]. Hence let S(t) = a{t)a{ty. We also assume controllability, 

that is given x,y G K", and f < h < t2 < T, 3u € ^2(^1,^2) (^-1) 

such that the solution Xt of Xt — A{t)xt + a{t)ut, satisfies Xt^ — x 

and xtr^ — y. We also assume that A{t),C{t),T,{t) are piecewise 

continuous, locally bounded functions of time t, and S(t) 't: for 

all t. 

Now consider the following optimal control problem. Let T < ti < T. 
We wish to maximize 

Jl'{xi,u)= lt{xt,ut)dt + (l)'{xT) (3) 

where Xt S ffi" satisfies 

Xt = ft{xt,ut) = A(t)xt + a{t)ut (4) 
xti = xi (5) 

and lt{xt,Ut) = ^xt C{t)xt — \\ut\^ and with u G L2{ti,T) and the terminal 
payoff is a quadratic in state variable, but parametrized by z G M" as a bivariate 
quadratic 

(j)^x) = (j){x, z) ^ ^x'Px + x'Sz + ^z'Qz (6) 

where P, Q are symmetric and S is invertible. 

The optimal control value function is defined to be 

V'{h,xi)^V^^{x^)= sup Jl\x,,u) (7) 

ueL2(tuT) 

for all X, z £ M" and ti G (T, T] . We shall show that above value function is a 
solution of a particular Hamilton- Jacobi-BcUman (HJB) PDE. 



Let W^{t,x) be the solution of following Hamilton- Jacobi-Bellman PDE on 
(T, T] X 

-^W'{t,x)^H{t,x,VW'{t,x)) (8) 

where VFF^(t, x) is the shorthand for -§^W^-{t, x), and with the boundary con- 
dition defined in (6) 

W''{T,x) = (l){x,z) VxeM" 

and 

H{t,x,p) = sup {p'ft{x,u)+lt{x,u)} 

= sup \p' {A{t)x + (T{t)u) + l^x'C{t)x 
= ^x'C(t)a; + x'A\t)p + ^p'S(t)j> 

where = cr(t)(T'(t). 

Lemma 2.1 For any z G R", there exists a solution to (8), (9), (10) in 
C°°{(f,T] X M"), and this is given by 

W^{t, x) = ^x'Ptx + x'Stz + ^z'Qtz 

where Pt, St, Qt satisfy Pt = P, St = S,Qt = Q and 

~Pt = A{t)'Pt + PtA{t) + C{t)+Pt^{t)Pt 

-St = {A{t) + ^{t)Ptyst (11) 

~Qt = St'nt)St 
and St is invertible for all t £ (T, T] . 

Existence of the solution Pt : — T < i < T is assumed in {AA). This combined 
with local boundedness, and piecewise continuity of coefficients guarantees the 
existence of St, and hence that of Qt for — T < t < T. Proof that it solves 
HJB PDE, is immediate by substitution in (8) and (9). Let us define, B{t) = 
-{A{t) + S(t)Pt). Then St^ = ^B{ti,T)ST, where $b is the state transition 
matrix of the system it — B(t)xt. By Abel-Jacobi-Liouville formula 

det^B{ti,T) ^ e-'*i ^' >0 

Since both $3(^1, T) and St — S are invertible, St^ = ^B{ti,T)ST is invertible 
as well, g 

Next we need a verification theorem to connect HJB PDE solution to the 
control value function. 



(9) 



^lul^dt 
2' ' 



(10) 



Theorem 2.2 Let x,z eM." and h e {T,T] . One has 

W''{ti,x)> Jl^"{x,u) yueL2{ti,T) 

and 

where Ut = u{t,xt) = a(tYVW^{t,Xt) — a{t)' {PtXt + Stz), which implies = 
and 

Vt{x, z) = Vt'{x) = W^-{x) = ix'PtX + x'Stz + ]^z'Qtz (12) 
Proof. Ijet u e L2{ti,T). 

Jl\x, u) = / {Lt{xt,ut) + {A{t)xt + CT(<)ut)'VM^^(t, a;*)) dt + (P{xt, z) 

T 

{A{t)xt + (T{t)ut)''^W''{t, xt) dt 

1 

which by definition of H 

<[ HixtyW'{t,xt)) dt + (j){xT,z)- [ {A{t)xt + <j{t)utyyW'{t,xt) 

which by (8) and (4) 

-^W^t, Xt) - it^W^t, dt + cf,{xT, z) 
^ d 



W''{t,xt)dt + (l){xT,z) 

dt 

=W'{ti,x) - WiT, xt) + 4>{xT, z) = W^^(ii, x) 

using (9). 

Also note that in the above proof, if we substitute ut = a{tyWW'^{t,xt), 
then using lt{x,u) = \x'C{t)x- \ut\'^/2, a{t)a'{t) = E{t) and (10), 

lt{xt,ut) + iA{t)xt + a{t)ut)'VW'it, Xt) 

= ^xt'C{t)xt + \vW^'{t, xt)ni)^W'{i, Xt) + A'{t)VW^{t, xt) 
= H{t,xt,VW'{t,Xt)) 

This converts the inequahty into equahty and we get J^^^{x,u) = W^{ti,x). ^ 

Now combining (7), lemma (2.1) and theorem 2.2 and substituting z — 
(zero vector), redefining Vt^{x) as Vix), and parametrizing the univariate value 
function by p instead of P for bivariate, we get the following corollary 



Corollary 2.3 Given t G (T, T] and x G M", the value function 
rT /I I,,. 12 



Vt{x)= sup \[ (\xtC{t)xt-'^-^^^^^-^ dt+]-XT'pxT)\ ^\x'ptx 

u<£L2(t,T) I Jti ^ / ^ J 

(13) 



subject to the dynamics 



xt = A(t)a;t + a{t)ut 

satisfies HJB PDF 

-§l'Ptix) = ^x'cit)x + x'A'it)vVtix) + ^vVt{xYj:it)vrtix) 

and pt satisfies the following DRE in which = (7(t)a{ty , 
-Pt = A{tYpt + ptA{t) + C{t) + pt^{t)pt 
with boundary condition px — P- 

This gives us the motivation to solve the DRE using underlying optimal control 
problem. The optimal control problem defined in (7) with bivariate quadratic 
terminal payoff parametrized by z will be useful in deriving the fundamental 
solution as will be covered in sections ahead. 



3 Fundamental Solution 

For given t2 and a general terminal payoff function (l){x) : M" M, let us define 
the operator, 

Sll\(i)]{x)= sup / lt{xt,ut)dt + (^{xt^) (14) 

ueL2{tid2) Jti 

Suppose that ti < t2 is such that the solution Sl^[(l)]{x) exists. (Is finite for any 
X € M".) 

We can restate (7) and (3) using above operator. Noting that Vrf{x) = 4>^{x), 
as defined in (6) , we have for all t € (T, T] 

v,^ix)^snr]{x)^snvs]{x) 

It is well known that operators S^^ form a semigroup. That is if ii < t < 
t2 < T, then 5j [0] = 5*J5*^[0]], which is precisely the celebrated Dynamic 
programming principle for this problem. That is with t2 = T, 

v;^(x) = simix) = si^ [sTm] (x) = si^ [F/] {x) 

/•* (15) 

= sup / lt{xt,ut)dt + Vt''{xt) 
ueL2{ti,t) Jti 



Maslov [21] proved that this semigroup is also Hnear in max-plus algebra. That 
is 

where (j)i{x), 4'2{x) are functions and fc G M is a scalar. If we define a® b = 
max(a, h) and a®h = a + b, then it is well known that (M U {— oo}, ©, (g)) forms a 
commutative semifield which is referred to as the max-plus algebra (see [18], [19], 
[20] for a fuller discussion) . 

We can extend it to functions so as to define the max-plus vector space. Let 
[a © b]{x) = max(a(a;), b{x)) and a{x) (g) A: = a{x) + k, where a, b : M" M. and 
A; S M. Using this notation, we have 

Si [01 ® 02] (x) = Si [0i] {x) © si [02] {x) 

Sl[k®(l)i]{x)^k®Sl[cj)i]{x) 

Now we shall define a max-plus kernel I : M" x M" M derived earlier in [2] 
and [8]. Let T < ti < t2 < T and x,y G M", and xt evolve with dynamics (4). 
Define 

Il2^^^y)^ h'''^ueut-^i.,y)lt^^t(^t,ut)dt iiUll{x,y)^% ^^^^ 
^ — oo otherwise 



where 



^tl{x,y) = {u e L2{ti,t2) : Xt, = x,Xt^ = y} 



Note that /jj^ — —oo indicates that it is impossible to reach y from x in time 
interval (ti,t2) using any possible control u. 

Fleming and McEneaney [2] proposed above kernel, and showed that 

Sl'^mx) = sup {Ill{x,y) + 0(y)) = / /*f (x,2/) ® ^(y) dy (18) 

j/eR" JR" 

and since depends only on the dynamics Xt — ft{xt, Ut) and running payoff 
lt{xt, Ut), it is independent of the terminal payoff ^(xtj). Hence it can serve as a 
Fundamental solution, and obtain Sl'^[<j)]{x) for any (f>{x) by a kernel operation. 

Remark 3.1 Note that Slf[(j)]{x) > -oo for all X £ M". Since for any u £ 
L2{ti,t2), solution to io = a; and Xt = ft{xt,ut) exists for t G [^1,^2] B.nd 
j^^ lt{xt,Ut) dt > —00. Hence 

Sll[4']{x) > lll{x,xtj + (f>{xt2) > / lt{xt,ut) dt + (j){xt2) > ~oo 

Remark 3.2 Also note that due to the controllability assumption {A.\), for 
ti < t2, we can always find control Ut which generates the trajectory x{t) satis- 
fying = X and Xt^ = y. Thus Ull (x, y) and ill v) - Iti h{xt,ut) dt > 
—00 for all x, y G M". For ti = ^2, lllix, y) — —00 for &\\y ^ x and /^^ (x, x) = 0. 



3.1 Computing the max-plus kernel 

First we derive a lemma about the end point of optimal trajectories. 

Lemma 3.3 Consider the system trajectory xt starting from xt-^ = x and evolv- 
ing according to (4) under optimal control Ut = a(ty{PtXt + Stz) from theorem 
2.2. Then for f <ti <t2 <T, 

St2S:t2 + Qt^z = St^'xt^ + Qt,z (19) 

Proof. By time- varying linear system theory, for a system evolving as per 

Xt = A{t)xt + cr{t)ut 

= A{t)S:t + cj{t)a{t)'{PtXt + Stz) 
= (A{t)+^{t)Pt)xt + ^(t)Stz 

solution is given as 

Xt2=^B{t2M)iu+ ^B{t2,T)Y.{T)SrZdT (20) 

Jti 

where <&s(t2,^i) = Ut^Uf^^, where Ut is the solution of differential equation 
C/t = B{t)Ut, with B{t) = A{t) + T.{t)Pt. 

It is well known that the state transition matrix 

*S(t)(i2,il)-*'_B(t)'(^l'^2) 

now, noting from (11) that St = -{A{t) + Y.{t)PtY St = -B{t)'St, and since St2 
is invertible, we have 

^B(t){t2,h) = ^'-BityihM) = {St,S^^^) ' = Si^^'St,' (21) 

Substituting in (20), and noting from (11) that ~Qt — St''S,{t)St, 

I't2 

*t2 = St^'Sti'xt^ + Si^' I Sr'^{T)SrZdT 

Jti 

= Si^^'St.'it, + ^t"'' (^^*' 5/S(T)5,dr^ z 
— ^t2^'^ti'xti + S^^^' {Qt^ — Qts) z 

thus we have, 

St2'xt2 + (9*2^ = St-^'xt-^ + QtiZ 



Remark 3.4 Note that Vz, since St^ and St, are invertible, (19) suggests a 
one-one and onto relation between start and end of optimal trajectories, xt-^ 
and Xt2. Thus Vy e M" there exists a x = St^' {St^'y + {Qti — Qts)^) such that 
optimal trajectory xt starting at Xt-^ — x, ends with y. Thus every y G M" is an 
optimal point for some initial condition. 



Remark 3.5 Note that due to max-plus linearity, if fc G M, using (15), 

V,\ = Si W + fc] (^) = si W\{x)^k = V,\{x) + k 

Thus while keeping the dynamics and the incremental payoff same, adding a 
constant to the terminal payoff only shifts the value function accordingly. The 
gradient hence the optimal feedback control remains the same. 

ut{x) = a{t)\/Vtl{x) - a{t)VVtl{x) = Ut{x) 

Hence the optimal trajectory, which is the solution to Xt = A(t)xt + a(t)ut{x), 
also stays the same. 

Now we shall prove another useful lemma before turning to the main result. 

Lemma 3.6 Given T < ti < t2 < T , and Qt evolving according to (11) with 
terminal value Qt — Q, then 

Proof. Note that since we assumed in {A.l) that the system if — A{t)xt + 
a{t)ut parametrized by {A{t), cr(i)) is controllable. This is true if and only if the 
following controllability grammian is invertible for any T < ti < t2 < T . 

[ ' ^A{ti,T)a{t)a{t)'<^A{tuT)' dt y (22) 

Thus for all x,y E M", 3 control Ut such that is the trajectory x — A{t)xt + 
a{t)ut with xt-i = X satisfies Xt2 — y- 

Now we claim that system {A{t) + I]{t)Pt, a{t)) is also controllable. This 
is clear because by using control ut — Ut — a{t)'PtXt, we can keep the system 
trajectory same and reach from x to y. 

X = A(t)xt + ij(t)ut 
^ {A{t) + a{t)a{ty Pt)xt + a{t) {ut - (T{t)' PtXt) 
= {A{t)+^{t)Pt)xt + <j{t)ut 

Hence similar to (22), using B{t) = A{t) + Y.{t)Pt and a{t)a{ty = E(t), following 
controllability grammian is invertible. 

[ \B{tur)i:{t)<^B{tur)' dt>Q (23) 

Substituting $s(ii,T) = S^_^^' Sr' from (21), 

f ' $B(i2,r)S(t)$B(t2,T)'dr = [ ' Si^^'Sr'nt)SrSy^^dT 
Jtl Jtl 

^ sr," (Qt, - Qt,) sr,' (24) 



where in last equation, we used Qt evolution from (11). Using (23) and since 
iSfj is invertible by Lemma (2.1), we have Qt^ — Qt2 ^ 0. ^ 



Theorem 3.7 Given x,y eM."- and T < ti < <2 < T, then 

inf [V,l (x) - V,l {y)] = III y) (25) 

Since by (12), Vi\{x) = \x' Pt^x + x' St^z + \z'Qt^z and V^l{x) = \x' Pt-^x + 
x' St^z + ^z'Qt^z, the max-plus kernel lll{x,y) is also bivariate quadratic. 

Ill{x,y) = ]^x' Iii\\x + x' Ii2\\y + ]^y' l22t\y where 

hit\ = Pti - St^iQti - Qti) ^St^' ^26) 
Ii2tl — Stj^iQti ~ Qt2) ^^t2' 

Proof. Let xt^ = x. Since E(t) > and St is invertible, by (11), Qt = 
-S't^{t)St ^ 0, hence Qt, - Qt^ h 0. For any z e M" 

Vt:{x)-Vtl{y) 

= Sll[V^Jx)-Vtl{y) 

= sup I / lt{xt,ut)dt + Vtl {xt2 ) - Vtl (y) 

substituting for V^^, 

f Z"*^ 1 1 

= sup II lt{xt,ut)dt+ -x't^Pt^xt^- -y'Pt2y+{xt2-yySt2z 

«ei2(ti,t2) Uti ^ ^ 

Since Ull{x,y) — {u <^ L2{ti,t2) ■ Xt, — x,xt2 = y} d -^^2(^1,^2) , and Vm S 
l^lli.x,y), Xt2 = y. 

f Z"*^ 1 1 

> sup II lt{xt,ut)dt+ -y'Pt2y- -y'Pt2y+iy-yySt2z 

= sup / lt{xt,ut)dt^ lll{x,y) (27) 

Taking infimum over all z G M", 

inf^ [t//^ (a;) - Vt\ {y)] > ill v) (28) 

Since Qt^ — Qt^ by 3.6, define z — {Qt, — Qt2)^^{St2y ~ St,'x). Hence 

St2'y + Qt2Z = St,'x + Qt,z 



hence using (19) the optimal trajectory xt starting from Xt^ — x and with 
terminal payoff V^^{-), ends at = y. Let the corresponding optimal control 
be Ut- Let us define k = —Vf^^{y) — — iJ^z'Qt^z + \y' Pt^y + y'St2z) to create a 
shifted terminal payoff function 

U,i{x) = V,l{x) + k=Vf^{x)-V,l{y) 

1 1 (29) 

= ^x'Pt^x - -y'Pt^y + ix- yySt^z 

From remark 3.5, Ut, Xt are also the optimal control and trajectory for the 
following problem with the terminal payoff . Hence 

Vf^ix) - Vl{y) = I sup /*' k{xuut) dt + V:^{x) \ ~ V^iy) 
yueL^itiM) Jti J 

= sup { r lt{xt,ut) dt + Vl {x) ~ Vl (y) 

M6L2(tl,t2) Wtl 

= sup / lt{xt,ut) dt + U^^ixt^) 
ueL2(tiM) Jti 

ltixt,ut)dt+ Ut^{xt2) 

ti 

since U^-^ixtJ = U^^{y) = from (29) and u G Ul^{x,y) 

< sup / lt{xt,ut) dt ^ lll{x,y) (30) 

Thus we have 

inf [V,\ (x) - V,l {y)] < Vl {x) - (y) < {x, y) (31) 

Hence (28) and (31) together give us (25) and also the following 

inf^ [V,l {x) - V,l (y)] = V,i (x) - Vl (y) = llf {x, y) (32) 

with z = {Qt, - Qt2)^^{St2'y - St/x). 

Substituting z in (32) and expanding, we get (26). □ 



Remark 3.8 It is interesting to note that the formulae extend graciously even 
when assumptions on controllability are violated. In that case Qti ~ Qt2 ^ 
and may not be invertible. We can do singular value decomposition 



Qu -Qt2^[ Ui U2 ] 



' A 


■ 




■ u[ - 













where [Ui U2] is unitary matrix, and A is diagonal matrix of nonzero eigenval- 
ues. We can obtain Moore-Penrose psuedoinverse as 

If we take limit of the formula (2G), by replacing all zero eigenvalue by fc > 0, 
and letting fc — > 0, then we obtain following formulae, which also give us a 
representation of the reachable set. 

(26) with {Qt, - Qt,y^ q 't ^ '1, P ranp-pCr/,^ 

7^*2/^ „\ _ J replaced by {Qt^ - Qt^) 



otherwise 



In the special case when E(t) is zero matrix, and no control is possible. Ui 
is empty, since there are no nonzero eigenvalues. (Qi — Q2) and (Qi — (52)^ 
are zero matrices. Hence range (J7i) ~ 0, which is the null range. Hence 
St^'x - St^'y = 0, and 

r x'Pt,x-y'Pt,y liy^Si^^'St.'x 
Ill{x.y)=l 

y —00 otherwise 

Thus using (21), only state accessible from starting point x is S^^'St^'x = 
^B{t){t2,ti)x = <^A{t){t2,ti)x, since B{t) = A{t) + S(t)Pt = A{t). This is the 
wellknown solution to time-varying linear differential equation, it ~ A(t)xt. 

Now we shall prove a theorem which can allow us to combine max-plus 
kernels in time. 

Theorem 3.9 Let T < ti < t2 < < T. then max-plus kernel llf can be 
computed from ll^ and I^^ as follows 

Illix.y) = Sll[lll{.,y)]{x) = sup {lll{x, z) + lll{z,y)] (33) 

zGK" 

Thus lll{x,y) = ^x' Iiillx + x' lutly + hv'^'^^tlv where 

^llti — ^llti - ^12ti (-12241 + ^11*2/ -'I2ti 
Il2lt^-Il2ll{l22ll+Illliy'll2ll (34) 
-(2241 — J22t2 - Jl2t2 (-12241 + JlltJ ^12*2 

Proof. Note that by remark 3.2, Z^*'(a;, y) ^ for all ta < tb and x,y e M". 
lll{x,y)^ sup / lt{xt,ut)dt 



since <^ (x, y) = U.^k^ «^ {x, z) n Ull (z, y)) 

rt3 



/■13 

= sup sup / lt{xt,ut)dt (35) 

Now, we consider the following 

sup / lt{xt,ut)dt 

u&ull(x,z)nUll(z,y) Jti 

< sup I / lt{xt,ut)dt+ / lt{xt,ut)dt 
ueu*^(x,z)nul^(z,v) wti "'t2 

< sup <^ / lt{xt,ut) dt> + sup <^ / lt{xt,ut)dt> 

ueUl'^{x,z) ^-^ti J ueUt'^{z,y) '-"'t2 J 

= /*,^(a;,z)+/*|(z,y) (36) 

Now, since ll^{x,z) > — cx), Ve > 0, 3 u S Ul'^{x,z) and trajectory xt with 
= a: such that 

Zt(xt,?2t)dt + e > /*^^(a;,z) (37) 

ti 

Similarly 3 {t G (z, y) and trajectory Xt with — ^ such that 

lt{it,ut)dt + t>lll{z,y) (38) 



t2 

Now we can create augmented control ii such that ttj = for t e [^1,^2) 
and wt = {t( for t G [^2,^3], and extended arbitrarily beyond. Note that if 
Xt is corresponding trajectory, then starting with xt^ — x, Xt = Xt for t G 
[^1,^2]- Hence Xt2 — z and Xt — Xt for t G [^2,^3], hence Xt^ — y. Hence 
u G Ull{x,z) f]Ull{z,y). Moreover using (37) and (38), 

/•ia r^s 

sup / lt{xt,ut) dt> / lt{xt,ut) dt + / lt{xt,ut)dt 



u£Ull{x,z)r\Ull(z,y)-'tt Jti Jt2 

lt{xt,Ut) dt + / lt{xt,Ut)dt 



tl Jt2 

= lll{x,z) + lll{z,y)~2e (39) 
Since e is arbitrary, from (36) and (39), we have 

/t3 
lt{xt,ut) dt = lll{x, z) + lll{z, y) 



which with (35) proves (33). Now, using (26) and since {Qt^ — Qt^) >~ and 

-'22ti + 

= -St, {{Qt, - Qt^r' + {Qt2 - Qt,r') St,' 
-< 

Thus {/*j*(a;, z) + lll{z^ y)} is concave in z. Thus supremum in (33) exists, and 
we get (34) by algebraic computation of the local maxima. □ 

Remark 3.10 Note that ll^{x^z) has the same bi-variate form as Vt given 
by (12), and both and Vt evolve in time interval (ii,i2) according to the 
semigroup 5*^ as per (33). Hence the parameters satisfy DREs similar to the 
(11). 

A(i)7n*^ + /iijM(t) + C{t) + /n*-'S(i)/n?^ 
{A{t)+mh,'^)'h2\' (40) 

/l2*^'E(t)/l2^ 

3.2 Algorithm 

Thus following is the final algorithm to obtain the fundamental solution, and to 
convert a particular solution of (11) into a general solution. It gives us closed 
form solution to the DRE (1) using max-plus kernel (18). We shall reiterate 
the formulae derived earlier to make the section self-contained. 

• Choose terminal t2 and the parameters {Pt,, St2,Qt2) of the terminal bi- 
variate payoff Vt,{x) = \x' Pt2X + x' St2Z + \z'Qt2Z, such that Pt2, Qt2 are 
n X n symmetric matrices, and St, is n x n invertible matrix. 

• Propagate (P, S, Q) backwards in time according to (11) till time ti < t2- 
That is 

-Pt - A{t)'Pt + PtA{t) + C{t) + Pt^{t)Pt 

~St^{A{t) + i:{t)Pt)'St 

-Qt = s/mst 

• Compute the max-plus kernel or the fundamental solution as per (26), 

Ilfix, y) = I^l{x, y) = \x'h,\\x + x' h2\\y + \y' h2\\y 



d 
dt 
d 

dt 

d_ 

~dt 



r t3 _ 



-'124 — 



^22t — 



parametrized by triplet {hi, I12, 122)11 "^t^^re 

hi — Pti — StiiQti — Qta) ^S**/ 

I12 = StiiQti — Qta) ^-Sta' 

h2 = ^Pt2 ^ St^iQti — Qti) ^StJ 

• As per corollary 2.3, if terminal payoff is given by Vt2{^) — \x'pt2X and 
if Pt evolves as per DRE (1), 

-p = A{t)'pt + ptA{t) + C{t) + pt^{t)pt 

and if exists, that is the solution does not blow up during ^2 ~^ ti 
evolution. Then by (18) 

Vt,{x) = Ix'pt.x = 5*f [nj(x) = sup Il'^{x,y)+Pt2{y) (41) 

Thus algebraically we get 

^x'pt.x ^ sup l^^x' Inx + x' Ii2y + ^y' {h2 + Pt2)y^ (42) 

Pt, = /ii - /i2(ft. + I22y^h2 (43) 

which is the analytical solution to the DRE (1). Thus we have converted a 
general solution to a bivariate DRE (P, S, Q) as per (11) into fundamental 
solution ill, ^'^'^ t^G"^ particular solution pt of (1). 

• As seen in (26), as (t2 — ti) 0, [Qt^ — Qti)^^ ™ay blow up, as Qti — > 
Qt2- Thus parameters of the max-plus kernel /^^ also blow up, causing 
numerical inaccuracy in propagation. To remedy this, and alternate form 
of propagation (43) is proposed as follows. 

After substituting kernel parameters from (26) in (43), with some manip- 
ulation we get 

= -(Qti - Qt2)~^ 

— [Qti — Qt2) ^ St2' {Pt2 ^ Pt2 ^ St2(Qti ~ Qt2)~^ St2') -^ta (Qtj ~ Qtg ) 

which by using Woodbury's matrix identity [22], 

= {~{QH-Qt2) + st2'{pt2- Pt2r^st2y^ 

Thus 

St,'{pt, - Pt^r'St, + Qt, = St2'{pt2 - Pt2r'St2 + Qt2 

or rearranging, propagation from to pt^ is given by 

Pti = Pt, - St, {Qt, - Qt2 - St2'{pt2 - Pt2r'St,y' St,' (44) 

This formula does not blow up for small time step propagation, and yields 
an accurate propagation. 



Remark 3.11 Note that we assumed that propagation ^x'pt^x = Slf[Vt2]{x) 
exists, and derived (43). This is also equivalent to I22 + Pt2 ~< 0; so that 
the supremum in (42) exists. Thus l22tl + ^2 ^0 characterizes all initial 
conditions pt^ for which DRE propagation t2 — > ti does not blow up. Also 
note that the minimum time T for which solution to DRE exists, depends on 
initial condition. Max-plus kernel obtained from one particular solution, may 
cause blow up for a different initial condition. Surprisingly it is possible to pass 
through singularity/ solution blow up using (43) without numerical instability 
(since instead of marching through singularity, we step over it), and generate 
solution trajectories akin to tan(a;) which is the solution to i: = 1 + x^. 

4 Semiconvex dual DRE 

Now we shall introduce the concept of semiconvex duality which can help us 
transform time invariant DREs into semiconvex dual DREs. 

4.1 Semiconvex duality 

A function p{x) : M" = MU {—00} is defined to be luniformly semiconvex 

with (symmetric) matrix constant K if V{x) + ^x'Kx is convex over M". We 
denote this space by . 

Semiconvex duality is parametrized by a bivariate quadratic kernel 

(f){x, z) — -X Px + x' Sz + -z'Qz (45) 

where P and Q are symmetric matrices. We use this kernel to define semiconvex 
duality. 

Theorem 4.1 Let V € S^^ , S is invertible and (t){x, z) defined as above. Then 
\/z e R" we can define the dual Q{z) of primal V{x) as follows. 

Q{z) = MiVix) - 0(a;, z)] = V4V]{z) (46) 

X 

from the dual Q{z), primal can he recovered again using 

V{x) = sup[0(x,z) + Q{z)] = V-^\Q]{x) (47) 

z 

(j){x,z) is called the kernel of duality. Thus 'D'^^'D ^\P\{x) — V{x). 
Proof. Since V e S^^ , 'P{x) — 4>{x, z) is convex in x. Now, 

sup[0(x, z) + Q{z)] = aup[4>{x, z) + iiii[V{y) 

z z y 

— sup inf [7^(2;) + (f){x,z) 

z V 

— sup inf [7^(2/) H — x' Px 

z V 2 



^ly'Py+{x-y)'Sz] 



Let z — Sz. Since S is invertible, z also spans K" 

= ]-x'Px + supinf [^'(y) - \y'Py + [x - y)'z\ 

= -x' Px + sup[a;'z + mi\P{y) y' Py — y'z] 

2 ^/ y 2 

Since 'P{y) — -^y' Py is convex, by Legendre-Fenchel transform (see e.g. Theorem 
11.1 in [27]) 

= -x'Px + Vix) - ^x'Px 
2 ^ ' 2 

= n^) 

□ 

If we choose V{x) — ^x'px, 4>{x,z) = ^x'Px + x'Sz + \z'Qz and assume 
p > P and S is nonsingular, then Vix) G S~^. Hence by substitution in (46), 
we get Q{z) = -^z'qz, where 

q^-S'ip-Py'S-Q (48) 
We can also derive following inverse relation 

p= -S{q + Qr' S' + P (49) 



Corollary 4.2 Using very similar methodology, if Q{z) + z'Qz is concave over 
z e W\ that isifq + Q^O, then 

V^V-'Q{z)^Q{z) (50) 

Remark 4.3 Let us observe that a result using (18) and (15) we saw earlier, 
can be reposed in the following manner using semiconvex dual notation. With 
ti <t2, since Sl^[(j)]{x) exists, and ll^{x,y) is bivariate quadratic, we have 

v,i {x) = sli [v,\]{x) = sup {III y) + K (y)) = ^^74 I^^^l (^) (si) 

y *i 

Since above supremum exists for all x, l22tl~^Pt2 ^ 0, hence by (43), and matrix 

congruence Pt, - hill = -^^^t {Pt, + h2tiy^ 112%' h 0. Hence V^'^ix) € 
5^^". Hence we can take semiconvex dual. Now using (50), 

inf {V,^^{x) - Illix, y)) = P,*. [V;^]{y) = V^.^V.] [V,l]{y) = V,l{y) (52) 

As a special case, from (41), we have 

VtAx) = SlllVt.Kx) = sup lll{x,y)+VtAy)^T^;.l[Vt,]{x) (53) 



4.2 Dual differential Riccati equation 

Now let us start with primal space quadratic V{x) = ^x'ptx, with pt varying 
with time as per (1). Let us find its dual Q{z) — ^z'qtz, with kernel 0(a;,z) = 
^x'Px + x'Sz + ^z'Qz. We assume pt — P> and S is nonsingular, so that we 
can use theorem 4.1. Using (48) and (49), we get 

qt^-S'{pt-Py' S-Q (54) 

Pt = -S{qt + Qr^ S' + P (55) 
Difi^erentiating both sides of (54), 

qt = s' {pt ~ py' Pt {pt -py^s (56) 

If the primal quadratic evolves according to (1), we can track the evolution of 
the dual. Substituting for pt from (1), and for pt from (55) in (56), we get 

qt Hqt + Q)S-' {A{typ+ PA{t) + PE(i)p) s^'\qt + Q) 

~{qt + Q)S-\A{t) + ^{t)Pys 
- S^{A{t) + ^{t)P)S^'\qt + Q) + S'^{t)S 
Using (11), and after simplification, we get. 

~qt - qtS-^PS-^^qt + qtS-\PS-'^Q - S) 

+{ps-^^Q - sys-^'qt + QS-^PS-^'Q 

-QS-^S~ {QS-^Sf + Q 
This show that the dual quadratic also satisties a Riccati equation 

-qt = A{tyqt + qtAit) + C{t) + qt'f^it)qt (57) 

with coeSicients 

A{t) = S-\PS-^^Q - S) 

S(t) = S-^PS-^'^ (58) 
C{t) = QS-^PS-^'Q 

-QS-^S-{QS-^Sf + Q 

where P, S and Q are constants defined by (11). Thus {A{t),f:{t),C{t)) = 
f{P,S,Q,A{t),m,C{t)). 

4.3 Kernel Matching conditions 

By using (97) and some algebraic manipulation, it can be easily shown that 
simultaneous equations (58) are equivalent to following simultaneous equations, 

^ A(t)'P + PA{t) + C{t) + PY.{t)P = SE{t)S' 

= {A{t) + S(t)P)'5 = S{-A{t) + S(i)Q) (59) 

= 5'S(i)S' = -A{tyQ -QA + C{t) + Qt{t)Q 



These give a neater feasibility condition for finding a kernel parameters, {P, S, Q) 
to transform Riccati equation (94) into any other Riccati equation (57). 

Remark 4.4 Observing the symmetry between the primal and dual DREs mo- 
tivates us to propose a dual problem with dynamics i = A(t)z + a(t)u, with 
aa{t)' — S(t) (if S > 0), and payoff lt{zt,ut) — ^z'C{t)z, and a corresponding 
dual semigroup 5*^ similar to (14). Using the symmetry of above equations, it 
can be easily proved that 

<j>t,{x,z) = Sll[cl>tA;m^) = ~Sll[-<l>t,{x,-)]{z) (60) 

Thus given the coefficients of primal and dual Riccati equations, both (Pq, Sq, Rq) 
and {Pt,St,Rt) satisfy (59), suggesting that these equations are not indepen- 
dent. Indeed, it can be verified that (59) are also equivalent to following matrix 
equation, which uses classic Hamiltonian matrices. 

ICn = niC (61) 



where 

K. is invertible, with 



s-^p -s-^ 

S'-QS-^P QS-^ 



A{t) E(t) 

-c{t) -A{ty 



-S + PS-^'Q PS-^' 



Hence H = JCTiJC ^ . Hence a necessary condition to find kernel (P, S", Q) to 
convert DRE {A{t) , C {t) ,i:{t)) into (^(t), C'(t), S(t)) is that matrices U and H 
be similar. Sufficiency conditions are being investigated. 

With K{t) = IC{P{t), S{t),Q{t)) and using (59), it is also easy to verify that 
-K: = K.H = UK. and = -JQ-^JCK.-^ = UK.-^ ^ JQ-^U Thus if $(ii,t2) 
is the state transition matrix associated with the linear time varying system 
i(t) = -n{t)x{t), then 

1Ct, - /Ct,$(i2,ti) and /C,-' = ^{t2,h)Kt^^ 

Similarly, if <&(ti,t2) is the state transition matrix associated with the linear 
time varying system i(t) — ~'H{t)x{t) , then 

= Ht2,ti)lCt, and /Cr,' = K^^^^h,h) 

If <I>(t2, ti) and $(^2, ^i) are partitioned into four n x n blocks, then we have 

Sr'Vu ^tT" 1 _ [ I'll '^12 1 [ Sr^^'Qt, 



-5*1 




■ $11 $12 ■ 








$21 ^"22 





-St2 



(62) 



and 



Matching terms, we get following set of equations 

St, = (*ii + <i>i2Pt,r"St, = St,{Qt,^i2 - 
Qt, = Qt, - 5t/($ii + ^i2Pt,r^'i>i2St, 

= (Qt2^12 + ^22y\Qt2^11 + ^21) 
Pt, = (*21 + *22PtJ (*11 + ^12Pt,V 
= Pt, - St,^i2{Qt,^12 + i'22)''^t.' 



^'22)^ 



fll ^12 
^•21 *22 

(63) 
(64) 
(65) 
(66) 



Remark 4.5 Note that (66) is equivalent to traditional Davison-Maki approach, 
substituting U = St^' andV = PtS^^', we have Pt = VU'^ and U and V 
evolve according to 







■ $11 


^'12 ' 






. y, . 




^'21 


^'22 




. y. . 



which is the solution to (2). 

We can substitute (64), (66) and (65) into (26), to get 

hltl = ($21 + *22PtJ(*ll + ^12Pt.,y^ + (*11 + ^12Pt,y^'^i2 
h2\l = -($11 + $12Pt,)'''$r2'(*ll + *12^tJ 
l22tl^~Pt,+<^l2i.'^ll+'^12Pt,) 

But since (/n, /12, ^22) depend only on {A{t), C(t), E(t), ii, ^2) and are indepen- 
dent of starting (P, S*, Q), we can take Pt^ = to simplify above formulae. 



/lit? = *2i*ir ' + *ir"*r2' 
h2\i = -<i>iri'<i>r2''i>ii 

l22\l = <i>r2''i>ii 



(67) 



Remark 4.6 Above formulae are useful in deriving analytical solutions for 
{P,S,Q). Especially for the time invariant case, <&(i2,^i) — e^^^'^^*!-* and 



Ht2,ti) 



^ p-n{t2-ti) 



But note that the eigenvalues of Hamiltonian 7i are 



symmetric along imaginary axis, thus contains both stable and unstable eigen- 
values. For time-invariant case, this leads to more and more ill conditioned 
4>(ti, 12), is thus useful as an analytic solution only for small t2 — ti. 



4.4 More Fundamental solutions 



Now we shall see how semiconvex duality can help us relate solutions of the 
primal and the dual DRE in various ways through max-plus kernel operations. 



in process generating other, possibly easier ways to compute fundamental solu- 
tions. 

With ti < <2 and u G ^2(^1,^2), first we define a backward trajectory of the 
system, with final point x, which is unique solution to 

xt = ft{xt,ut) t e [^1,^2], with a;t2 = y (68) 

which is guaranteed by assumptions {A.V). 

Theorem 4.7 With T < ti < t2 < T , backward dynamic programming coun- 
terpart of (15) also hold true. That is if Xt is the backwards trajectory ending 
at Xt2 — y, under controls u, as defined in ( 68), then 

Vtliy)^ inf {- T' k{xuUt)dt + V,\{xt,)\ (69) 
= ini^{V,\{x)-lll{x,y)) (70) 
^Sll[V,\]{y) 

Note that this also defines the semigroup operation Sl^ for ti < t2. Also we 
have V,l = Sll [V,\] - S*^^ [Sll [F/J] . Thus 

Sll=Sll-^ (71) 

Proof. Given u E ^2(^1,^2), let Xt be the backwards trajectory which satisfies 
(68) with Xt2 — y. Let x — Xtj^- Then x ^ y — J^^ ft{xr, Ur) dr, thus 

/'t2 rt 
Xt=y- / ft{Xr,Ur)dT = X + / ft{Xr,Ur)dT 



Thus Xf is also unique solution of xt = ft{xt, Ut) for t E [ii, t2] with Xt-^ — x. 

Thus we have y — Xt2, where it — ft{xt,ut), with xt-^ = x. Hence using 
(15), 

Vtl{x)> r It{xt,ut)dt+V^^{xt,)= f'lt{xt,ut)dt + Vtliy) 
Since x = Xt^, 

Vtl{y)<Vtl{xt,)- rk{xt,ut)dt 
taking infimum Vm G ^2(^1,^2) and corresponding backward trajectories, 



< inf - / k{xu Ut) dt + Vt\ {xt,) (72) 

ueL2{ti,t2) I Jt^ 



Now specifically, we can take x — S^^' {St^'y + {Qt^ — Qu)) and the forward 
trajectory xt starting from xt-^ — x as per optimal feedback control discussed 
in lemma 3.3, Ut = u{t)'{PtXt + Stz). By remark 3.4, it is clear that Xt^ = y. 
Hence 

V^\{x)= sup / It{xt,ut)dt + Vtl{xt2) 

lt{S:t,Ut)dt+Vtl{xt,)^ / lt{S:t,ut)dt + Vtl{y) 



Since Xt-^ = x, 



Vtliy)^- r lt{S:uU,)dt + Vl{xt,) 



> inf -/ ltixt,ut)dt + Vt'^{xtJ\ (73) 
and we prove (69) using (72) and (73). We get (70) from (69) and (52). □ 



For the forthcoming analysis, we assume following. 

Assume that ti < t2- Vtix) = ^x'pt{x), withpt evolving as per (1). 

We also define a bivariate quadratic function (^^{x) = 4>tix,z) = 

^x'PtX + x' Stz + ^z'Qtz, with parameters {Pt, St, Qt) evolving as 

per (11). A duality kernel (/>^(x) = 4>l^{x) = (fit^ix, z) = \x' Pt^x + (^.2) 

x'St2Z+ ^z'Qt^z transforms Vt into Vip^Vt = Qt = -^x'qtix). Also 

assume that pt2 h Pt2-< St2 nonsingular, and that Vti{x) = 

Sll[Vt2\{x) and cpl{x) = Sll[Vt2]{x) exists. 

Theorem 4.8 Assuming {A.2), the semiconvex dual of Vti{x) under kernel 
(j)t-^ (x, z) exists, and is same as the semiconvex dual of Vt2 (x) under kernel 
(j)t2{x,z). Using (46), Vz e K" 

inf [Vt2 {x) - (bt2 {x, z)] = inf [Ft, (y) - <f>t, (y, z)] = Q^, (z) (74) 

X y 

That is 'Dff,^^ \Pt2] — ^0ti ["^til- terms of parameters, following equation holds 
true. 

qt2 = -St2'{Pt2 - Pt2)~^St2 - Qt2 = -St/iPti - Pti)~^Sti - Qti (75) 

Proof. Existence of the dual I?^^ [T'tj] is evident from theorem 4.1. Note 
that from corollary 2.3 and (18), with x — xq, 

\x'pt^x = VtAxti) = [T't^KxtJ sup {l*tl{x,y) +Vt2{y)) 



using (43), = /ii -/i2(pt2 +I22) ^h2 ■ Since Vt^ exists, pt^ +I22 is invertible. 
Similarly since exists, = 15^*2 . Thus 



2 ^ 

exists, thus by a similar logic. In — Pt^ is also invertible. 
Now we shall consider a bivariable function 

^{x, y) = -ci>i (x) + III y) + (y) (76) 

note that given x G M", 



argmax y) = argmax if^ (x, y) + T'tj (y) 
V y 



(77) 



= (/22 + Pt2)"'^i2'a: 
similarly given y E M", 

argmax '0(a;, y) — argmax —0^^ (x) + ij^ {x, y) 

X X 

= a.rguimlllix,y)-(j>l{x) (78) 

X 

= -{lii-Pu)-'li2y 

Using 3.6 and (26), Qti—Qt2 >~ and hence /12 >- 0. Combined with invertibility 
of /22 + Pt2 and III — Pti- from (77), (78), we have 

y = {y ■ y — argmax -(/^(x, y) for some x E M"} = M" 
J/ 

A" = {a; : x = argmax -0(2;, y) for some j/ e M"} = M" 

a: 

Hence by corollary A. 2 from appendix, 

inf sup ip{x,y) — inf sup ip{x,y) (79) 
y y X 

With this preparation, we are ready to prove the main result, 
inf [Pt, (x) - (x)] = inf [Sll [Vt^Mx) ~ cj^^ {x,, )] 

= inf [sup {III V) + {y)) - ^t, {x) 

X \_ y 

= inf sup {^cbi {x) + III y) + (y)) 

X y 

= inf sup ip(x,y) 

X y 



using (79), 



inf sup xp{x,y) 
y X 



= inf 

y 



7't,(2/)+sup(-0^^(x)+/,\^(a;,y)) 

X 

VtAy)--^nf{ct>l{x)-ilf{x,y)) 



inf 

y 



Algebraically, it is easy to see that, minimum occurs at a; = {pt-^ — Pj J~^S't^z 
and y = {pt^ — Pt2)~^St^z, respectively. Plugging this into (74) gives us (75). 
Note that these equations have the same form as (44) obtained earlier, g 

Corollary 4.9 Assume (^-2). Using (f>(x,z) — 4>t^{x) as a duality kernel, let 
\z'qtz — Qt{z) — V^Vtix) for all t g [ii,t2]- The semiconvex primal of Qt2{z) 
under kernel 4>ti {x, z) is same as the semiconvex primal of Qt^ {x) under kernel 
(t>tAx,z). Using (46), \/x G M" 

sup [QtM + K z)] = sup [Qt, (z) + {x,z)\= Vt, (80) 

z z 

That IS V-} [Qt,] = V-} [Qt,], and 

Pt, = -StA<it2 + Qt.r'St/ + Pt, = -St^iqt, + Qt^r^Stj + (81) 

Proof. We have, 

Q,^{z) =V^[Vt,]{z) ^V^^^[Vt,]{z) (82) 
Qt2{z) = V4Vt2]{z) = V^jVt^Kz) (83) 

Also using theorem 4.8, V^^^ [Pt2]{z) = 'D^^^ [PtiK^), hence 

J^tJ] = ^ti (84) 

Thus using, (82), (83) and (84) 

%l [QtJ = V, [^^'2 [^*^]] = = V. [^^'. [^*^]] = [2*=] 
Finally we obtain (81), using (80) and (47). □ 



Now we shall obtain on a time-varying version of the result previously ob- 
tained in [1], to complete our picture of kernel relationships between primal and 
dual DREs. For this result, we make an additional assumption. 

We assume that F{Pt) = A{t)' Pt + PtA{t) + C{t) + PJl{t)Pt = 

~Pt >- for t e [^1,^2]- Thus we have, pt^ - pt^ ^ 0. ^ ' ' 



Qt2 {z) under kernel -B^J (z, z) 



Theorem 4.10 Assuming {A.2) and (A. 3), Qti{z) is the semiconvex primal of 
QtA^)= sup [Bll{z,y) + QtM] or 

QtA^) = 'D-'[Qt,]{z) 



(85) 



where, 

^tl (z, y) = inf {0*1 {x, y) - (j)t^ (x, z)} (86) 

Hence Bll{z,y) = \z' B^^\\z + z' B^2ty + h' B22ty , with 

Bull ~ -^St^iPti — Pti) ^St2 — Qt2 
B,2ll^St2'{Pt,-Pt2)-'St, (87) 

and 

qt, = Bn^: - Bi2^:(B22?; + qtJ-'Bi2ll' (88) 
Proof. In corollary 4.9, we saw that V'^^ [QtJ — T^^^ [Qtal- hence we have 

Q^Az)^V^,V^l[Qt2\{z) 

= inf feMQtJ(x)-0i,(a;,z)) 



inf sup ((/)ti (x, y) + Qf^ (y)) - (a;, z) 
inf sup ((/>t, (x, y) + Qt., (y) - (I>t2 {x, z)) 



— inf sup ^p{x,y) (89) 



where 

V'(a;, y) = ^2:'(Pfi - Ffja; + x'St.y + ^y'{Qt2 + ^2)2/ - x'St^z - ^z'Qt^z 

Note that by (A. 3), Pt^ — Pt^ >- 0. Hence ip{x,y) is strictly convex in x. Also 
observe that by corollary 4.9, Vt^x) — supy{(f>tj^{x,y) + Qt^iy)) exists for any 
X e M". Hence Qtj + qt^ -< 0. Thus il}{x,y) is strictly concave in y. For such 
a convex-concave function following saddle point exists. By setting V and 
Vj,?/; equal to zero, and solving, we get 

= ((F,, - Pt^) - 5i,(gt, + gtj-^^t/)"' 5t,z 

For such x^ and y*^, ip{x^,y) < '0(a;'^, ?/'') < ip{x,y^). Hence by a well known 
result, 

inf sup ip{x,y) — ^{x'^ ,y'^) = sup inf 'tp{x,y) (90) 



Using (89) and (90) 

Qti(z)= sup inf ip{x,y) (91) 

= sup (2*2(2;)+ inf {(l)t,{x,y) - (j>t^{x,z))] (92) 
= sup {QtM + Bll{z,y)) (93) 

(87) can be easily obtained from (86) by finding local minimum in x (which 
is global minimum, since infimum exists), substituting and term- wise equating 
coeSicients. Similarly (88) results from substituting Qt = ^z'qtz , (86), (87) 
into (85). □ 



Thus equations (43), (53), (75), (74), (81), (80), (85), (88) can be summa- 
rized in the diagram below. Note that primal and dual quadratics are on top 
and bottom respectively. Vertical and diagonal lines show duality transforma- 
tion with indicated kernel. Arrows are directed from the primal to it semiconvex 
dual. 

Thus in conclusion, so far we saw three distinct ways of solving (1), that is 
obtaining pt^ from pt^ ■ 

1. Direct method which assumes only (A.l). Formulae are given by (43) and 
(26). Propagation is achieved by following transform. 

Problem with this method is that as — > ^2, parameters of the kernel 
blow up, limiting solution accuracy. 

2. Alternate method, which assumes (A.l) and (^.2). Formulae are given 
by (75), which is same as (44). Propagation is achieved by following 
transform. 

This method works better for small time step propagation, since parame- 
ters of kernels and (f>t2 do not blow up. 

3. Third method assumes (^.1), (^.2) and (^.3). Time invariant version 
of this method was first proposed in [1]. Formulae are (87) and (88). 
Propagation is achieved by following transform 




Problem with this method is similar to the direct method. Namely, ti 
t2, parameters of the kernel I^.^ blow up, limiting solution accuracy. 



5 Time Invariant problem 

The theory so far developed for DREs with time- varying coefficients, extends 
readily for time invariant DRE, in which A{t) = A, C{t) = C, a{t) = a hence 
S(t) = aa' = S. Again we assume (AA). We state such DRE again for 
reference. 

-p = A'p + p + C + p'Y.p (94) 

If we define (x) = ^x'pt^x, then with T < ti < t2 < T, l{xt,ut) = ^Xt'Cxt — 
^|utp, and dynamics Xt = f{xt,ut) = Axt + crut and starting point Xt^ = x, 
value function is 



Vt,{x) ^ Sl![Vt,]{x) ^ sup I l{xt,ut)dt + Vt,{xtj] (95) 

Then Vt^{x) = }^x'pt^x, where pt satisfies DRE (94). Let us define /S. — t2 — ti. 
Let (5 G M. Note that using the time invariance of dynamics and incremental 
payoff, and change of variables t ^ {t ~ 5), 

Sll[(j)]{x) ^ i sup / l{xt,ut)dt + (j){xt2)] ■■ xti = x\ 

= < sup / l{xt,ut)dt^<^{xt2^i)\:xt^-f,^x\ 

[«eL2(ti-<5,t2-<5) j J 

= sizimx) 

As a special case, using 8 = — i\ and 8 ~ t\ — ti respectively, 

0*2 cO CO C*2~*l — 

'-'tl — '-'ti-t2 ^ '-'-A — '-'O — Oq — O 

By similar argument, the fundamental solution, or max-plus kernel 

III = ^tVt2 = ^-A = i'r'^ = It = (96) 

Similarly, the bivariate quadratic DRE in (11) turns into 
-Pt - A'Pt + PtA + C + Pt^Pt 

-St = {A + ^Pt)'St (97) 
—Qt — St'^St 



5.1 A doubling algorithm 

Using (34), we can derive a useful doubling algorithm for solving (94). First we 
propagate the triad {Pq, Sq, Qq) backwards in time by A to obtain (P_Ai ^-a, Q-a), 
using (97). Then we can comput as defined in (96). Thus using (26) with 



ti = -A and t2 = 0. 

^11^ = -P-A — <S'-a(Q-a - Qo)^^S-a' 

/i2^ = 5_A(g-A - Qor'S„' (98) 

122'^ = —Pq — '5'o(Q-A — Qo) ^So' 

We can build up using (33) as follows 

I'^ix, y) = I^^ix, y) = sup {l^{x, z) + /i^(z, y)} = sup {l^{x, z) + /o^(z, y)) 

Z Z 

Thus using (34) 

r2A_TA rA/rA,T A\-1t A/ 
Jll — Jll — J12 (-(22 +^11 ) -'12 

^2^'^ = ^Il2^ {122'^ + ^11^) ^12^ (^^) 

r 2A _ T A r A//T A, r A\-lr A 

J22 — -'22 ~ J12 (^22 + Jll j J12 

Thus by recursively combining kernel operations, we can get ^ m. N steps. 
Hence the name doubling algorithm. Finally the DRE (94) can be solved by 
(43). This can give exponentially fast propagation of the DRE. 

Remark 5.1 In [1], an analogous doubling algorithm was derived in semiconvex 
dual space. We shall sketch it here without proof. It starts with construction 
of the dual kernel i3*J = using time invariant special case of (87). 

Sii^ = -5o'(P-A - Por^So - Qo 

Bi2^ = So'{P-A - Por'S^A (100) 

^22^ = -S-a'{P-A - Por^S-A + Q-A 

and kernel is doubled using following formulae. 

r) 2A _ p; A o A/o A.r) A\-l p Ai 
-Dll — -Dll — i'12 (i'22 + iSll ) -012 

_Bi2^^ = — i3l2^ (-622^^ + -Bll^) -812"^ (101) 

p 2A_p A o Ai frj A.r) A\^l n A 
-D22 — -D22 — i'12 (£>22 + -Dll ) ^12 

Thus for T = 2^ A, is found using N doubling operations, and propagation 
is achieved using V-t = 'D^g'^BT'D^f,[Po]- Which using (46) and (47) implies 



P-T = -So [B^i - Bl, {Bg - So'ipo - Por'So - Qo) Pf^' + Qoj So'+Po 

(102) 

Remark 5.2 For both (98) and (101), as A ^ 0, (Q_a - Qo)^^ and (P_a - 
Po)^^ can be grow very large, leading to numerical inaccuracy at very small 
time step propagation. A work-around is designed in the next subsection. 



5.2 Alternate doubling algorithm 

First we deduce following generalization of theorem 4.8. 

Lemma 5.3 Consider t < 0. Define (pf{x) = (f>t{x,z) — ^x' PfX + x' Stz + 
^z'QtZ, with parameters {Pt, St, Qt) evolving as per (97). Assume that a primal 
function Vq{x) is such that Vq{x) — x'PqX is a convex function. Also assume 
that Vt{x) = St\PQ\{x). Hence a duality kernel (f>^-{x) = 0o(^) ~ 0o(2;, z) = 
^x' Pqx + x' Sqz+^z'Qqz transforms Vt into T^^^Vt = Qt- Then, the semiconvex 
dual ofVt under kernel 4'q{x, z) exists, and is same as the semiconvex dualVo^x) 
under kernel 4>q{x, z). 

inf [Vo{x) - ^)] - inf [Vt{y) - Mv, z)] = Qo{z) (103) 

X y 

Proof. The proof proceeds exactly as in theorem 4.8 with t2 = and ti = t. 
The only difference in assumption is, instead of assuming Vtix) — ^x'ptx with 
Po '>~ Pqi now we assume 7^0(2;) — x' P^x is a strictly convex function, and Vt = 
StVo. This enables broader choice of Vq, which we shall use in next corollary. 

□ 

Corollary 5.4 Letti < andt2 < 0. Staring with terminal (t>o{x, z) we assume 
that 4>ti{x, z), (^(2(0;, z) and (pti+t^ix, z) exist, and 4>t2{x, z) — (j)Q{x, z) is strictly 
convex. With choice ofVo{x) = (^^^(a;) = (x.z), and t = ti, we have 

Vt = Sl[4>l]{x) = 5,\VtJ0y(a;) = <PUt2i^) = K+t^i^.-A 
Hence using lemma 5.3 for all z and z 

inf (x, z) - (t)o{x,z)] = mi[(j)t^+t2{y,z) - (j)t^{y,z)] (104) 

X y 

Hence following is true. 

Pti+t2 = Ptl + Sti {{Qo — Qti) + So'{Pt2 — Po) ^Sq} St/ 

St,+t2 - St, {(Qo - QtJ + So'{Pt2 - Po)-'^o}"' So'{Qt2 - Qo)-^St2 (105) 

Qti+t2 — Qt2 + St2' {{Po — Pt2) + SoiQti ^ Qo) ^Sq'} St2 

Proof. Substituting the parameters of cj)t for t = 0,ti,t2,ti + ^2 m (104), 
using (48) and (49) and matching terms gives us following 

So {Pt2 ^ Po) ^Sf) + Qo = Sti {Pti+t2 ^ Ptl) ^Sti + Qti 

St2'{Pt2 ~ Po) ^Sq — St-^+t^ {Pti+t2 ~ Ptl) ^St-, 
— St2'{Pt2 ~ Po) ^St2 + Qo = ^Sti+t2 {Pti+t2 ~ Ptl) ^Sti+t2 + Qti+t2 

from which (105) can be derived by some manipulation, and using Woodbury's 
matrix inversion formula. 



Corollary 5.5 Using ^1=^2 = —A, A > 0. We have following kernel doubling 
formula, in terms of parameters of the biquadratic duality kernel. 

P-2A = P-A + S-A {{Qo - Q-a) + So'iP-A - Por'Soy' 
^-2A = 5_A {(Qo - Q-a) + So'{P-A - Poy'Soy' So'{Q-A - Qo)~'^-A 
Q-2A = Q-A + S-a' {{Po - P-a) + So{Q-A - Qor'So'y' 5_A 

(106) 

Using these we can create kernel (f>t{x,z) parametrized by {Pt, SitQi) at time 
t = —2^ A in N steps and achieve Riccati propagation using time-invariant 
version of (44), that is 

Pt = Pt- St {Qt - Qo - So'ipo - Por'Soy' St' (107) 



5.3 Dual DRE and Analytic solutions 

Extending (57) to time-invariant case, we have if Qt{z) = V^lVt] = \z'qtZ, 
where (j){x, z) = ^x'PtX + x'Stz + ^z'QtZ, then qt satisfies following dual DRE. 



It 



A'qt + qtA + C+qt'^qt 



(108) 



Similarly extending (59) to the time invariant case, coefficients of dual DRE, 
{A, C , E) satisfy following matching or compatibility conditions. 



= A'P + PA + C + PT.P = 5SS" 



-S^{A + Y.Pys = S{-A- 
-Q = S'Y.S^-A'Q-QA~ 

Similarly extending (61), 

/C7i — TiK, 



C + QtQ 



(109) 



(110) 



where 



s-^p -s-^ 

S'-QS-^P QS-^ 



A 

-c 



-A' 



A_ 
-C 



-A' 



Note that with a constant duality kernel, hence constant /C (assuming invert- 

ibility), and using time invariance oi H, H — IC'HIC~^ = 0. Hence dual DRE is 
also time invariant. 

Extending the figure 5.3 to the time varying case gives us figure 4.4. Note 
that ti = —A and t2 = 0, and (f)t{x,z) = ^x'PtX + x'Stz + \z'Qtz, with the 
triad {Pt,St, Qt) evolving as per (109). Note that we have one new relationship 
in figure (4.4), stated as below. 



Lemma 5.6 



2?0a[^o] = Q-A 



Proof. From figure (5.3), we have [PtJ = V~l V^^^ [Vt^ = Sll [Pf J. Because 
of time invariance, for any (5 G M, 

tl <Pti ti+d <Pt2+S <Pti+S 

Hence using S — A, ti — —A and t2 — 0, we have, 

V-A = Sl^iVo] = S^[Po] - V^^V^jVo] 
Thus finally, using P^JT'-a] = Q-a, 

V^^Vq = V4,„[V-a] - Q-A 

□ 

Remark 5.7 Note that letting A oo gives us relation between the solutions 
of primal and dual Riccati equations. That is, if p is stabilizing solution of 
= A'pt+PtA+C+pt'Spt and q is the stabilizing solution of A'qt+qt-A+C+qtf^qt. 
Then q is the semiconvex dual of p, under kernel (P, S,Q)o- Hence 

q = -So\p-Po)-'So-Qo (111) 

Remark 5.8 Note that -^z'qtz — Qt — 'D^g\Pt\{z) and qt satisfies dual DRE. 
Similarly, using lemma 5.6, 'Dif,^\Pt\{z) = Qt-A = \z'qt-AZ. Since qt~A also 
satisfies same dual DRE (108), kenel 0a (a;, z) also satisfies same matching con- 
ditions (109). Let us define 

St' — QtSf Pt QtSi 

Thus from figure (4.4), if /Co satisfies matching conditions (110), then so does 
/Ca, for any A > 0. 

The kernel matching conditions (59) allow us to transform one Riccati equa- 
tion into any other. In particular, we can transform a Riccati equation into a 
Linear equation amenabe to analytical solution. Thus we can derive well known 
analytical solutions to Riccati equations derived earlier in [4] and [6]. 

5.3.1 Analytical solution 1 

We wish to solve pt = A^pt + ptA+ C + pt^pt- To convert this into a linear 
equation, let us choose Sq = I , C = — Q. Thus kernel matching conditions 
(59) imply that, Vt 

-Pt ^ A'Pt + PtA + C + Pt^Pt^Q (112) 
-St = {A + ^Pt)'St = -StA (113) 
-Qt = S'^^St^-^Qt-QtA (114) 



From (112), for all t, Pt is the solution to the Riccati equation A^P + PA + 
C + P'SP — 0. Thus if the Riccati equation has stable and unstable solutions, 
say P- and F+, and if Pq = P+, then Pt = P+ Vt. 

From (113), since So = I, A = -{A + T.Pf. Solving Sq ^ I and S ^ -SA, 
we have S-t = e~^*. 

From (114), Qq is solution to the lyapunov equation, A'^Qq + QqA + E = 0. 
It can be proved that Qq = (P+ — P_)^^. Solving Q = —A^Q — QA, we get 
Q.t = e-^"*goe-^*. 

Combining all above and (44), we get 

{p-t ~ P-t) 
- -S-t [-5'o'(po - PoY^Sq + Q-t - Qq] 

-Atr I r, \-l I -A^t^ -At 1-1 -K^t 

Rearranging and substituting Qq = (P+ — P-)^^ and A = —(A + EP^)-^, we 
get 

= {{VQ-P^)-^-{P--P^)-^\ (115) 

Equation (115) is same as the Method 3 developed by Leipnik in [4]. Note 
that (115) assumes that extremal solutions P_|_ and P_ are well separated. 

This method being analytic, is extremely fast and accurate, and does away 
with the need to combine operators. 

5.3.2 Analytical solution 2 

Above method suffers from numerical errors when extremal solutions of the 
Riccati equations, are not well separated. Following almost analytic method 
can solve this problem. The reason it is called almost analytics is because it 
requires computation of integrals of matrix exponentials instead of exponentials 
themselves. 

If we use Sq = I, = 0, Qo = 0, then matching conditions (59) give us 

-Pt ^ A'Pt + PtA + C + Pt^Pt^O (116) 
-St = iA + J:PtySt = -StA (117) 
-Qt = S'J:St = -A'Qt-QtA + C (118) 

Hence again, Pq ~ Pt = P for all t, where P is any solution of Riccati equation 
(116). S-t = e^*-^, with -A = (A+ SP)"^ as earlier. Since Qq = 0, C = E 
and Qt evolves according to (118). If T is a solution of lyapunov equation, 
-A'T + TA=C, then Q^t = e-*^'Te-*^ - T. Hence 

-Qt = -A'e-'^'Te-'^ -e-'^'Te-'^A 

= -e-*^'(A'r + r^)e-*^ 

-tA' /=i -tA 
— e Ce 



Since Qo = 0, and C = E, Q_t = /J e~*'^'Ee"*'^. Substituting in (44), and 
using -A = {A + SP)^ = B, 



JB' 



JO 



dr 



(119) 



This formula works well for problems with extremal solutions not well sep- 
arated. In limiting case, when P = P^ = P_, it can be proved that —A = 
{A + SP+)^ = 0, which imphes 

(-p+p_0"' = (-^ + po)-'-a 

Using shorthand for pt ~ pt — P and pQ ^ pq — P It can be shown that above 
is equivalent to 



„tB' 



(120) 



This formula is same as one developed by Rusnak [6] . Integrals of matrix expo- 
nentials can be calculated easily using techniques from [7]. For our use, if 



exp 



B 



Flit) F2{t) 
F^{t) 



then /J e*^ te^* dt = Fsity F2{t) . 



5.4 Numerical Experiments and Results 

Thus we have seen three different approaches to solve the time invariant DRE, 
(94) in this text. Following is the summary of algorithms. 

• Starting with {Po, So,Qo), find the solution of bivariate DRE (97) on a 
smaller time interval, [— 1,0]. This can either be done analytically using 
(64), (65), (66) (for a small time interval, since Davison-Maki theorem is 
illconditioned for large t), or by time marching using appropriate solver. 
Here we use fixed step Runge-Kutta fourth order method for the same. 
We shall denote number of steps by Nrk. Note that an often useful initial 
condition is Pq — 0, Sq = I and Qo — 0. 

• Using {P-t, S-t, Q-t), we can construct kernels P* as per (100) and /* as 
per (98). We may also choose to find / analytically using (67). 

• Having found a kernel at time —t and starting with given po, we can choose 
to perform AI kernel doubling operations to solve for kernel at time 
and N time-stepping operations to compute p^T for T ~ 2^^ Nt. We can 
do this three different ways as follows. 

Method A Construct * using doubling formula (99) recursively M 
times, and evolve po back in time to get p^T using N stepping oper- 
ation as per (43). 



Method B Construct * using doubling formula (101) recursively M 
times, and evolve po back in time to get p_T using TV stepping oper- 
ation as per (102). This method was first proposed in [1]. 

Method C Use doubling formula (106) M times recursively to compute 
(P, Q, S)_2Mf. Evolve po back in time to get p-T using N stepping 
operations as per (107). 

Assuming that we evolve the solution from to t, using Runge-Kutta fourth 
order method, using Nrk steps, and perform M doubling and N stepping opera- 
tions, computational complexity (flops) needed for three methods are as follows 
(found using [26]) 



Method A Nrk{32n^+3n^)+{M+l) (l6n2 + f _^ ^n)+N {iSn^ + fn^ + ^ 
Method B Nrk{32n^+3n'^)+{M+l) {l6n'^ + + ^n)+N (54n2 + 19n^ + 11 
Method C Nrk{32n^+3n^)+M {Gln^ + f -I- fn)+N {36n^ + f -|- fn) 



Now we shall apply the theory developed so far to the example problems. 
We shall also benchmark Method B proposed in [1] against Method C, which 
show greater accuracy at very small time step propagation. 

5.5 A stiff time invariant example 

We shall show in this section that algorithms discussed so far are applicable 
to stiff DREs. Usual time marching methods are constrained in step size by 
the stability requirement. To get a rough idea,let us look at a linear system 
x{t) = Ax{t). Explicit methods like Euler method impose condition |l-|-/iA| < 1 
on stepsize h, where A is an eigenvalue of A. If A is real and negative, this 
condition implies < j^. Thus excessively small step sizes may be needed to 
ensure accuracy if A is large. Following numerical experiments show that new 
algorithms give accurate answers with step size h = significantly larger 
than 2/ A when the initial transient phase is almost over. 
We shall choose example 4.1 from [11]. Symmetric DRE 



where diag[i] denotes the diagonal matrix with successive diagonal entries 1, 2, . . . , n 
and U is an orthogonal matrix. Above DRE is solved for different n (size of the 
equation) and k (larger the value of k, the stiffer the equation). The error of 
the computed solution is defined as err = \ \p — p\\ p / \\p\\ p , where X denotes the 
true and X denotes the computed solution. The analytical solution of above 
DRE is 



p{t) = ~p^{t) + ein; 



p{0) = Udiag[i\U 



t > 



p{t) = [/diag 



fcsinh kt + icosh kt 



U 



-1 



cosh kt + ^sinh kt 



where diag[-] here denotes the diagonal matrix whose diagonal entries are gen- 
erated by letting i take successively the values 1,2, ... ,n. 



Let < ti < t2- Starting with p-ti = p{^ti) which is the analytic solution, 
we computed the solution p-t2, for different ti,t2,n, k, and using various choices 
for Nrk, M, N. Note that solution is marched for t = T/{N2'^), and results 
tabulated. 



tl 


tl 


n 


k 


M 


N 


Nrk 


h 


2 

|A| 


err 





10 


10 


5 


1 


10 


10 


1 


0.2 


2.94 X 10-^ 





1000 


10 


5 


10 


1 


10 


1000 


0.2 


1.50 X IQ-^ 





1000 


10 


10 


10 


1 


10 


1000 


0.2 


8.07 X 10-6 





1000 


100 


10 


10 


1 


10 


1000 


0.02 


8.07 X 10-6 





1000 


100 


10 


10 


1 


100 


1000 


0.02 


1.55 X 10-12 





1000 


100 


10 


10 


1 


500 


1000 


0.02 


3.24 X 10-15 





1000 


100 


10 


16 


1 


10 


1000 


0.02 


6.65 X 10-1° 





100000 


100 


10 


22 


1 


_* 


100000 


0.02 


2.98 X 10-15 





100000 


100 


100 


22 


1 


_* 


100000 


0.02 


2.24 X 10-" 


3 


10 


100 


100 


8 


3 


10 


2 


0.02 


1.52 X 10-6 


3 


10 


100 


100 


8 


3 


1 


2 


0.02 


1.77 X 10-2 





10 


100 


100 


11 


1 


100 


10 


0.02 


9.84 X 10-12 





10 


100 


100 


1 


2" 


100 


10 


0.02 


4.83 X 10-9 



* Note that a blank entry for Nrk implies that the solution for {P, Q, S)t was 
found analytically using (64), (65), (66). 

Usually error in the final solution arises out of 



1. Error in initial time marched solution (P,S,Q)t and hence in /*, which 
gets amplified due to doubling algorithms. Higher steps in time marching, 
Nrk, lead to a more accurate final solution. 

2. Error in computing /* for very small t. Since matrix inverses involved in 
computation of /* in (99), grow more singular as t ^ 0. For extremely 
small timesteps, parameters of /* can blow up, leading to inaccuracies. 
Hence it is preferable to keep t above a certain minimum threshold. 

Exact error analysis of these methods is beyong scope of this paper. But above 
results show that method 1, is fast, useful and numerically stable fundamen- 
tal solution for long time horizon propagation of a stiff DRE. Kernel doubling 
methods give a quadratic convergence rate, rather than the linear rate of step- 
ping iterations. This enables us to take longer time steps, much higher than the 
stability threshold of the stiff DRE explicit time marching algorithms. 

5.5.1 Time invariant doubling algorithms 

Now we shall bechmark the numerical performance of method 1, 2 and 3 using 
the example from [1]. We consider a 2 x 2 size DRE, with 

0.216 -0.008 " 
-0.008 0.216 



A = 



-2 1.6 
-1.6 -0.4 



,C = 



1.5 
0.2 



0.2 

-0.4 



starting with initial condition 



Po = 



-0.1 






-0.1 



and we use the duality kernel {P, S, Q)o = [D, —D, D), where 

D = 



-1 -0.2 
-0.2 -0.1 



We also compute a near exact solution using variable step Runge-Kutta 
(ode45 in MATLAB) method for T = 4, with absolute tolerance 10~^^ and 
relative tolerance of 10"^'^. Treating this as the true solution, we compare 
this against solution by all three methods, using just one step of Runge Kutta 
fourth order fixed step method for computing {P, S,Q)-t and using only kernel 
doubling. Thus we choose TV = 1, Nrk = 1 and vary M between 3 and 17. Error 
in the computed solution is defined as err — \\p — p\\ p / \\p\\ p ■ Here p indicates 
the true solution, and p is the computed solution. Figure 3 show plots of M 
against err for all three methods. We can see that the method 3, can achieve 
maximum accuracy of 9.48 x 10^^'^ as against to 5.77 x 10^^° for method 1 and 
3.85 X 10^^" for method 2 first proposed in [1], which is significantly better. 
Although increasing M , decreases err initially, eventually all method exhibit an 
increase in err with M. This is because t — 4/2*^ 0, and — Pq 0, 
S-t — 5*0 — > and Q-t — Qo — > 0, and their inverses grow very large and can 
not be accurately computed. Thus it is advisable to keep t above a minimum 
threshold. In this example, with M = 13, i = 4.88 x lO^''. 

We also benchmark above methods against a fourth order Runge-Kutta 
method, for various stepsizes, and plot err against computational cost mea- 
sured in flops (using [26]) in figure 4. Specifically, the flopcount for Runge 
Kutta fixed step 4th order method was taken to be Nrk{32n^ + 3n^). Note that 
above methods exhibit a faster convergence rate, and yield more accuracy for a 
given computational cost. 



6 Conclusion 

Thus in this paper, we derived a new fundamental solution for the DRE (1) 
with the time varying coefficients. We showed its equivalence with the max- 
plus based fundamental solution (17) and (18), proposed in [2], which arises out 
of a two-point boundary value problem and dynamic programming. We proved 
that such a fundamental solution is a bivariate quadratic for the linear/quadratic 
problem in (26). 

We showed that under a semiconvex duality kernel, a primal DRE is trans- 
formed into a dual DRE. We derived the compatibility conditions which the 
coefficients of such DREs should satisfy in (59). These conditions can be read- 
ily expressed as er (61) in terms of classic Hamiltonian matrices of system dy- 
namics and symplectic matrices of kernel parameters. This led us to derive 



analytical forms of the kernel parameters involving state transition matrix of 
the Hamiltonian system as per (64), (65), (66). 

Using backward dynamic programming principle, we showed the kernel rela- 
tionships between time shifted primal and the dual DREs, shown in Figure 4.4. 
This enabled us to develop three different methods to propagate a DRE. Time 
invariant analogue of one of these methods was first proposed in [1]. 

As a special case of the time-varying problem, we derived doubling and 
stepping algorithms in primal space, for the time-invariant DRE ((98), (99), 



Using the kernel matching conditions, we could design kernels which con- 
vert primal DREs into linear differential equations under duality, which have 
wellknown analytical solution. Using this, in sections 5.3.1 and 5.3.2, we re- 
discovered the analytic solutions for the DRE, previously developed by [4] and 



Finally in section 5.4, we used these methods on the stiff problems and 
demonstrated their speed, accuracy and numerical stability. We also bench- 
marked them against previous dual space doubling algorithm developed in [1], 
and demonstrated an order of magnitude improvement in the accuracy for prop- 
agation at a very small time step. 

Thus this paper provides an elegant fundamental solution for the time- 
varying DRE, useful for stiff problems and long time horizon propagation. It 
also provides a powerful unifying framework based on optimal control formula- 
tion, semiconvex duality and max-plus algebra, which enables us to solve Riccati 
differential equation, and see existing methods in new light. 

I would gratefully like to acknowledge Mr. Ali Oran for introducing me 
to this problem, and my advisor. Prof. William McEneaney for his inspiration, 
encouragement and fruitful discussions. I would also like to thank Prof. William 
Helton for his guidance and help on NCAlgebra, which was useful in discovering 
some of the formulae described here. 

Appendix 

Lemma .1 Let (j){x,y) he a bi-variable function of x lE M" and y G M'" . Let us 
define sets 



(43)). 



[6]. 



y — {y' : y' — argmax(/)(a;, y) for some x G M"} 

V 

X — \x' : x' = argmax0(a;, y) for some y G M™} 



(A-1) 



(A-2) 



X 



Then following inequalities hold true 



inf sup (j){x,y) < inf sup 4>{x,y) 

y&y xeK.-^ 2;eR" j^gRm 



(A-3) 



inf sup 4'{x,y) < inf sup 4>{x,y) 



(A-4) 



Proof. Note that (A-4) follows from (A-3) using symmetry in x and y. 
To prove (A-3), take any yi G y. From definition of 3xi E M" such that 
yi = argmaxj^gjj™ (^(x, y). Hence 

sup > = max y) 

> inf sup <p{x,y) 

Taking infimum over all yi E y, we get (A-1). Hence proved. ^ 

Corollary .2 If (f){x,y) is a bi-variable function of x E M" onrf y G I^"- M^i^ft 
A" and y as defined in (A-1) and (A-1) respectively, if X = y = M", then using 
(A-3) and (A-4), we have 

inf sup <j){x, y) = inf sup (f){x, y) (A-5) 
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Primal DRE: -pt = A{typt + ptA{t) + C{t) + pjl{t)pt. 




Dual DRE: -qt = A{t)'qt + qtA{t) + C{t) + qtf:{t)qt. 
Figure 1: Time varying problem: Duality relationships. 

Primal DRE: -pt = A'pt +ptA + C + pt^pt- 




Dual DRE: -qt = A'qt + qtA + C + gfEgt. 
Figure 2: Time invariant problem: Duality relationships. 




Figure 3: N vs. Solution Error 




Figure 4: N vs. Solution Error 



